plus sn10x, SI

Baf53cre ENS neurons, SI data
CR infection 7d, CTL x4 CKO(Ahr-cko) x4

loading 60k cells,
cellranger called 35,059 cells

pure neurons downstream

load dependancies

load obj

GEX.seur <- readRDS("./GEX0112.seur.pre.rds")
GEX.seur
## An object of class Seurat 
## 24081 features across 18878 samples within 1 assay 
## Active assay: RNA (24081 features, 1084 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap
color.FB <- ggsci::pal_d3("category20c")(20)[c(1,6,11,16,
                                               2,7,12,17
                                              
                                               )]

color.cnt <- color.FB[c(1,5)]
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "preAnno2") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

filtering

GEX.seur <- subset(GEX.seur, subset = preAnno1 %in% levels(GEX.seur$preAnno1)[1:4] & DoubletFinder0.1 == "Singlet")
GEX.seur
## An object of class Seurat 
## 24081 features across 16562 samples within 1 assay 
## Active assay: RNA (24081 features, 1084 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "preAnno2") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

markers.new.s <- list(EMN=c("Chat","Bnc2","Tox","Ptprt",
                       "Gfra2","Oprk1","Adamtsl1", 
                       "Fbxw15","Fbxw24","Chrna7",
                       "Satb1","Cntnap5b",
                       "Gabrb1","Nxph1","Lama2","Lrrc7",
                       "Ryr3","Eda","Tac1",
                       "Kctd8","Ntrk2","Penk",
                       "Fut9","Nfatc1","Egfr","Mgll",
                       "Chrm3"
                       ),
                 IMN=c("Nos1","Kcnab1",
                       "Gfra1","Etv1",
                       "Man1a","Airn",
                       "Adcy2","Cmah","Creb5","Vip","Pde1a",
                       "Ebf1","Gpc5"
                       ),
                 IN=c("Npas3","Synpr","St18","Gal",
                      "Kcnk13",
                      "Neurod6","Moxd1","Sctr",
                      "Piezo1","Sst","Adamts9","Kcnn2",
                      "Calb2"),
                 IPAN=c("Calcb","Nmu","Adgrg6","Pcdh10",
                        "Ngfr","Galr1","Il7","Aff2",
                        "Gpr149","Cdh6","Cdh8",
                        "Clstn2","Ano2","Ntrk3",
                        "Cpne4","Vwc2l","Cdh9","Scgn",
                        "Vcan","Cck","Piezo2","Kcnh7",
                        "Rerg","Bmpr1b","Skap1","Ntng1",
                        "Tafa2","Nxph2"),
                 CONTAM=c("Actl6b","Phox2b"),
                 pDEGs=c("Grid2","Ide","Btaf1","Cdh10","Slc15a2","Ccser1","Hdac9","Col25a1","Rspo2"))

pm.All_new.s <- DotPlot(GEX.seur, features = as.vector(unlist(markers.new.s)), group.by = "preAnno2",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All")
pm.All_new.s

pm.CTL_new.s <- DotPlot(subset(GEX.seur,subset=cnt %in% c("CTL")), features = as.vector(unlist(markers.new.s)), group.by = "preAnno2",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="CTL only")
pm.CTL_new.s

pm.CKO_new.s <- DotPlot(subset(GEX.seur,subset=cnt %in% c("CKO")), features = as.vector(unlist(markers.new.s)), group.by = "preAnno2",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="CKO only")
pm.CKO_new.s

re-clustering

GEX.seur <- FindVariableFeatures(GEX.seur, selection.method = "vst", nfeatures = 1500)

# Identify the 10 most highly variable genes
#top10 <- head(VariableFeatures(GEX.seur), 10)
top20 <- head(VariableFeatures(GEX.seur), 20)

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(GEX.seur)
plot2 <- LabelPoints(plot = plot1, points = top20, repel = T)
plot1 + plot2

head(VariableFeatures(GEX.seur), 200)
##   [1] "Mgat4c"        "Zfp804a"       "Cntn5"         "Gal"          
##   [5] "Adgrg6"        "Gm30613"       "Kcnb2"         "Klhl1"        
##   [9] "Gm32647"       "Robo2"         "Cntnap2"       "Cdh8"         
##  [13] "Sst"           "Cpne4"         "Ntng1"         "Nrxn3"        
##  [17] "Gm29521"       "Gm21847"       "Ano2"          "Pdzrn4"       
##  [21] "Brinp3"        "Lingo2"        "Cdh9"          "Ntrk3"        
##  [25] "Gpr149"        "Ebf1"          "Nmu"           "Kcnip4"       
##  [29] "Nrg1"          "Sgcz"          "Gm15261"       "Tmeff2"       
##  [33] "Astn2"         "Dgkg"          "Zfp804b"       "Pcdh9"        
##  [37] "Cdh18"         "Cmah"          "Prkg1"         "Clstn2"       
##  [41] "Car10"         "Schip1"        "March1"        "Cdh6"         
##  [45] "Grm5"          "Gpc5"          "Nxph2"         "Vip"          
##  [49] "Piezo2"        "Chsy3"         "Sema5a"        "Fgf14"        
##  [53] "Kcnq5"         "Asic2"         "Vwc2l"         "Rasgef1b"     
##  [57] "Plxna4"        "Egfem1"        "Dcc"           "Grin3a"       
##  [61] "4930509J09Rik" "Unc5d"         "Pbx3"          "Ccbe1"        
##  [65] "Efr3a"         "Pcdh10"        "4930432L08Rik" "Itgb6"        
##  [69] "Cysltr2"       "Zfhx3"         "Spock3"        "Cpa6"         
##  [73] "Gm20754"       "Gm15680"       "9530059O14Rik" "Col25a1"      
##  [77] "Csmd3"         "Myl1"          "Penk"          "Pde1a"        
##  [81] "Kcnh7"         "Robo1"         "Dgkb"          "Tnr"          
##  [85] "Serpini1"      "Trhde"         "Cacna2d3"      "Gm26612"      
##  [89] "Dpp10"         "Acp7"          "Skap1"         "Wdr17"        
##  [93] "Fut9"          "Plcl1"         "Luzp2"         "Gm15584"      
##  [97] "Csmd1"         "Abca9"         "4930422I22Rik" "Nxph1"        
## [101] "A330008L17Rik" "Gm49938"       "9530026P05Rik" "1700034P13Rik"
## [105] "Pde4d"         "Tafa2"         "Gm11342"       "Gm20713"      
## [109] "Bmpr1b"        "Tac1"          "Cntn4"         "Il1rapl1"     
## [113] "Trps1"         "Gna14"         "Pcdh11x"       "5530401A14Rik"
## [117] "Cux2"          "Bmp5"          "Cdh19"         "Arhgap6"      
## [121] "Gm29771"       "Kcnd2"         "8030451A03Rik" "Galnt18"      
## [125] "Chrm3"         "Kif26b"        "C730002L08Rik" "Dgki"         
## [129] "Gm16685"       "Kctd16"        "Prune2"        "D030068K23Rik"
## [133] "Kctd8"         "Hs6st3"        "B230323A14Rik" "Rab27b"       
## [137] "Myh3"          "mt-Co3"        "Calcb"         "Ghr"          
## [141] "Slc4a4"        "1700111E14Rik" "Lama2"         "Grik1"        
## [145] "mt-Atp6"       "Gm49127"       "B230312C02Rik" "1700063D05Rik"
## [149] "Mme"           "Iqgap2"        "Ano5"          "Galntl6"      
## [153] "Nos1"          "Gm11339"       "Ttc29"         "4930587E11Rik"
## [157] "Chrna9"        "Olfr1369-ps1"  "Gm12239"       "Gm20560"      
## [161] "Gm26632"       "mt-Co1"        "Gm16226"       "Galr1"        
## [165] "AC150683.1"    "Gm30382"       "Rasgrf1"       "Gm28494"      
## [169] "Fbxl7"         "Gm49678"       "Gm40841"       "Rxfp1"        
## [173] "C79798"        "Rerg"          "Meis2"         "mt-Co2"       
## [177] "Gm32916"       "Grm7"          "1700003I16Rik" "Adam2"        
## [181] "Cadm2"         "Adamts9"       "Synpr"         "Flrt2"        
## [185] "Alcam"         "Thsd7b"        "Cntn3"         "Nell1"        
## [189] "Met"           "Thsd4"         "Plxdc2"        "Tbx20"        
## [193] "Sorcs3"        "4931415C17Rik" "Dock10"        "Cdh10"        
## [197] "Ephb1"         "Dapk2"         "4833423E24Rik" "Mid1"
GEX.seur <- ScaleData(GEX.seur, features = rownames(GEX.seur))
## Centering and scaling data matrix

PCA

# exclude MT genes  and more 

#DIG <- grep("^Tra|^Trb|^Trg|^Trd|^Tcr|^Igm|^Igh|^Igk|^Igl|Jchain|^Hsp|^Rps|^Rpl|Hbb-|Hba-|^Dnaj|^AY|^Gm|^Hist",rownames(GEX.seur),value = T)

DIG <- grep("^mt-|^Tra|^Trb|^Trg|^Trd|^Tcr|^Igm|^Igh|^Igk|^Igl|Jchain|Mzb1|Vpreb|Lars2|Jun|Fos|^Hsp|^Rps|^Rpl|Hbb-|Hba-|^Dnaj|^AY|^AC|^AA|^AW|^BC|^Gm|^Hist|Rik$|-ps",
            rownames(GEX.seur),value = T)
CC_gene <- Hmisc::capitalize(tolower(as.vector(unlist(cc.genes.updated.2019))))


VariableFeatures(GEX.seur) <- setdiff(VariableFeatures(object = GEX.seur),
                                                c(#MT_gene,
                                                  DIG, 
                                                  CC_gene) )

GEX.seur <- RunPCA(GEX.seur, features = VariableFeatures(GEX.seur))
## PC_ 1 
## Positive:  Ntrk3, Ano2, Tmeff2, Robo2, Cdh8, Cpne4, Nrxn3, Mgat4c, Adgrg6, Clstn2 
##     Cntn5, Pdzrn4, Myl1, Gpr149, Spock3, Zfp804a, Cdh6, Cysltr2, Pcdh10, Cacna2d3 
##     Grin3a, Sgcz, Dgkg, Astn2, Iqgap2, Kif26b, Plxna4, Arhgap6, Cux2, Itgb6 
## Negative:  Lrrtm4, Galntl6, Tox, Ndst4, Grik1, Bnc2, Pcdh7, Cacnb2, Tshz2, Cdc14a 
##     Kcnd2, Epha5, Chrna7, Specc1, Lrp1b, Galnt17, Syt6, Plcxd3, Zfp536, Fbxw15 
##     Kcnab1, Synpo2, Brinp2, Nos1, St6galnac3, Lrrc4c, Rspo2, Fbxw24, Ptprt, Man2a1 
## PC_ 2 
## Positive:  Rbfox1, Bnc2, Ptprt, Tafa1, Gpc6, Tshz2, Grik1, Frmd4b, St6galnac3, Fbxw15 
##     Brinp2, Tcf7l2, Cdc14a, Plcxd3, Tox, Fbxw24, Tmem132c, Dock2, Pcdh7, Pld5 
##     Caln1, Oprk1, Unc5c, Specc1, Adamtsl1, Slc35f4, Stxbp5l, Nfia, Ccbe1, Sdk2 
## Negative:  Nos1, Egfem1, Auts2, Kcnq5, Cadps2, Etv1, Gfra1, Epha5, Dach1, Asic2 
##     Kcnab1, Alcam, Schip1, Dgkb, Cmah, Kcnt2, Stxbp6, Rgs6, Hs6st3, Tmem108 
##     Adgrl3, Lrrc4c, Creb5, Cdh11, Fat3, Il1rapl1, Kcnj3, Ebf1, Tenm3, Cacnb2 
## PC_ 3 
## Positive:  Sgcd, Adgrg6, Cysltr2, Gpr149, Nos1, Nmu, Slc4a4, Grin3a, Nfia, Dapk2 
##     Itgb6, Fgf13, Ano2, Ccbe1, Rgs6, Dgkg, Cbln2, Kcnab1, Lhfpl2, Ngfr 
##     Cpne4, Otof, Zfp536, Syt2, Gfra1, Efr3a, Zfp804a, Pex7, Calcb, Smad6 
## Negative:  Cdh18, Kcnip4, Csmd3, Klhl1, Kctd8, Cadm2, Pbx3, Gpc6, Htr4, Cntn3 
##     Meis1, Dlc1, Gabrg3, Skap1, C79798, Csmd1, Tenm2, March1, Car10, Serpini1 
##     Pde10a, Sema5a, Vwc2l, Khdrbs2, Dmd, Pde4d, Cdh9, Sybu, Sv2c, Pakap 
## PC_ 4 
## Positive:  March1, Klhl1, Sema5a, Vwc2l, Cdh9, Pbx3, Kcnh7, Prune2, Zfhx3, Chsy3 
##     C79798, Ntng1, Rasgef1b, Tenm4, Galnt18, Bcl2, Dcc, Galnt17, Zbbx, Ebf1 
##     Il1rapl1, Lncbate1, Olfr78, Kcnb2, Alk, Sez6l, Scgn, Bmpr1b, Kif26b, Etv1 
## Negative:  Dock10, Lingo2, Kcnip4, Ndst4, Nxph1, Tac1, Gda, Sorcs1, Prkg1, Epha5 
##     Csmd3, Cntn5, Thsd4, Dmd, Lrrtm4, Lrrc7, Kctd8, Kcnt2, Lrp1b, Fgf13 
##     Cntn3, Lama2, Ccdc60, Nyap2, Hgf, Pld5, Nhs, Ryr3, Calcrl, Grem2 
## PC_ 5 
## Positive:  Trhde, Cntn4, Ebf1, Trps1, Nrg1, Ptprd, Cpa6, Csmd1, Zmat4, Col18a1 
##     Kcnd2, Rmst, Egfem1, Gal, Chsy3, Myo1e, Nkain3, Npas3, Luzp2, Nav2 
##     Shisa6, Fstl4, Moxd1, Sctr, Prkg2, Myo16, Ntng1, Adgrl2, Asic2, Prkd1 
## Negative:  Dgkb, Kcnt2, Thsd7b, Klhl1, Alcam, Prkg1, Vwc2l, Rasgef1b, Il1rapl1, Fgf13 
##     Pbx3, Lrrc4c, Cdh20, Cdh9, Dpp10, Vcan, Khdrbs2, Zfhx3, Mgat4c, Galnt18 
##     Stard13, C79798, P3h2, Auts2, Zbbx, Slc44a5, Olfr78, Nmur2, Slc2a13, Rgs6
length(setdiff(VariableFeatures(object = GEX.seur),
                                                c(#MT_gene,
                                                  DIG,CC_gene)))
## [1] 1046
setdiff(VariableFeatures(object = GEX.seur),
                                                c(#MT_gene,
                                                  DIG,CC_gene))[1:300]
##   [1] "Mgat4c"   "Zfp804a"  "Cntn5"    "Gal"      "Adgrg6"   "Kcnb2"   
##   [7] "Klhl1"    "Robo2"    "Cntnap2"  "Cdh8"     "Sst"      "Cpne4"   
##  [13] "Ntng1"    "Nrxn3"    "Ano2"     "Pdzrn4"   "Brinp3"   "Lingo2"  
##  [19] "Cdh9"     "Ntrk3"    "Gpr149"   "Ebf1"     "Nmu"      "Kcnip4"  
##  [25] "Nrg1"     "Sgcz"     "Tmeff2"   "Astn2"    "Dgkg"     "Zfp804b" 
##  [31] "Pcdh9"    "Cdh18"    "Cmah"     "Prkg1"    "Clstn2"   "Car10"   
##  [37] "Schip1"   "March1"   "Cdh6"     "Grm5"     "Gpc5"     "Nxph2"   
##  [43] "Vip"      "Piezo2"   "Chsy3"    "Sema5a"   "Fgf14"    "Kcnq5"   
##  [49] "Asic2"    "Vwc2l"    "Rasgef1b" "Plxna4"   "Egfem1"   "Dcc"     
##  [55] "Grin3a"   "Unc5d"    "Pbx3"     "Ccbe1"    "Efr3a"    "Pcdh10"  
##  [61] "Itgb6"    "Cysltr2"  "Zfhx3"    "Spock3"   "Cpa6"     "Col25a1" 
##  [67] "Csmd3"    "Myl1"     "Penk"     "Pde1a"    "Kcnh7"    "Robo1"   
##  [73] "Dgkb"     "Tnr"      "Serpini1" "Trhde"    "Cacna2d3" "Dpp10"   
##  [79] "Acp7"     "Skap1"    "Wdr17"    "Fut9"     "Plcl1"    "Luzp2"   
##  [85] "Csmd1"    "Abca9"    "Nxph1"    "Pde4d"    "Tafa2"    "Bmpr1b"  
##  [91] "Tac1"     "Cntn4"    "Il1rapl1" "Trps1"    "Gna14"    "Pcdh11x" 
##  [97] "Cux2"     "Bmp5"     "Cdh19"    "Arhgap6"  "Kcnd2"    "Galnt18" 
## [103] "Chrm3"    "Kif26b"   "Dgki"     "Kctd16"   "Prune2"   "Kctd8"   
## [109] "Hs6st3"   "Rab27b"   "Myh3"     "Calcb"    "Ghr"      "Slc4a4"  
## [115] "Lama2"    "Grik1"    "Mme"      "Iqgap2"   "Ano5"     "Galntl6" 
## [121] "Nos1"     "Ttc29"    "Chrna9"   "Galr1"    "Rasgrf1"  "Fbxl7"   
## [127] "Rxfp1"    "C79798"   "Rerg"     "Meis2"    "Grm7"     "Adam2"   
## [133] "Cadm2"    "Adamts9"  "Synpr"    "Flrt2"    "Alcam"    "Thsd7b"  
## [139] "Cntn3"    "Nell1"    "Met"      "Thsd4"    "Plxdc2"   "Tbx20"   
## [145] "Sorcs3"   "Dock10"   "Cdh10"    "Ephb1"    "Dapk2"    "Mid1"    
## [151] "Gabrg3"   "Cbln2"    "Tenm2"    "Hs3st4"   "Npas3"    "Galnt13" 
## [157] "Ppp3ca"   "Stk31"    "Mast4"    "Loxl1"    "Arhgap15" "Cadps2"  
## [163] "Olfr78"   "Ror1"     "Gpm6a"    "Prr16"    "Dach1"    "Rbpms"   
## [169] "Tenm4"    "Rmst"     "Tenm1"    "Zfp286os" "Zbbx"     "Hccs"    
## [175] "Ngfr"     "Aff2"     "Nox3"     "Ntm"      "Calcrl"   "Wbp2nl"  
## [181] "Etv1"     "Sctr"     "Mecom"    "Hcn1"     "Col8a1"   "Pla2r1"  
## [187] "Tafa1"    "Rarb"     "Otof"     "Adamts6"  "Nell1os"  "Rtl4"    
## [193] "Sorcs1"   "Atg4a"    "Cpne8"    "Lrrc4c"   "Dab1"     "Vcan"    
## [199] "Wee2"     "Kcna1"    "Nwd1"     "Auts2"    "Serpine2" "F13a1"   
## [205] "Stxbp6"   "Kirrel3"  "Bfsp2"    "Necab1"   "Abca8a"   "Itgbl1"  
## [211] "Fstl4"    "Lhfpl2"   "Nek1"     "Ptprt"    "Pantr1"   "Il18r1"  
## [217] "Btk"      "Grp"      "Slc44a5"  "Slc2a13"  "Serpinb7" "Serpinb5"
## [223] "Nkain3"   "Gabrb1"   "Xlr4a"    "Ptk7"     "Opcml"    "Rbfox1"  
## [229] "Moxd1"    "L3mbtl4"  "Igsf21"   "Rgs6"     "Prkag2"   "Fpr1"    
## [235] "Col18a1"  "Gfra1"    "Lmcd1"    "St18"     "Brinp2"   "Slc44a4" 
## [241] "Edil3"    "Naaladl2" "Slco1a6"  "Ropn1"    "Slc6a7"   "Bcl2"    
## [247] "Scgn"     "Casp8"    "Ptprd"    "Grid1"    "Lrp1b"    "Htr4"    
## [253] "Cntnap5b" "Nfatc1"   "Tex15"    "Aff3"     "Bnc2"     "Scn7a"   
## [259] "Shisa9"   "Gda"      "Tacr1"    "Inhbc"    "Prkca"    "Gpc6"    
## [265] "Esr1"     "Pakap"    "Runx1"    "Kcnn3"    "Hs3st2"   "Ntrk2"   
## [271] "Apol7e"   "Slco1a1"  "Zbtb7c"   "Siah3"    "Mgat4e"   "Hsd17b3" 
## [277] "Lncbate1" "Eya1"     "Bche"     "Spag16"   "Plcb1"    "Adgrl2"  
## [283] "Lrrc43"   "Vmn2r65"  "Tmprss4"  "Antxr2"   "Pde10a"   "Scnn1a"  
## [289] "Scube1"   "Lin9"     "Pde7b"    "Nsun2"    "Htr3b"    "Apol7b"  
## [295] "Gata4"    "Olfr92"   "Sv2c"     "Adcy2"    "Cask"     "Dmd"
DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "orig.ident") +
  DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "orig.ident")

DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "FB.info", cols = color.FB) +
  DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "FB.info", cols = color.FB)

DimHeatmap(GEX.seur, dims = 1:12, cells = 1500, balanced = TRUE,ncol = 4)

decide PCs to use
ElbowPlot(GEX.seur,ndims = 50)

PCs <- 1:18
GEX.seur <- FindNeighbors(GEX.seur, dims = PCs, k.param = 20)
## Computing nearest neighbor graph
## Computing SNN
GEX.seur <- FindClusters(GEX.seur, method = 'igraph' ,resolution = 2)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 16562
## Number of edges: 641121
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8276
## Number of communities: 26
## Elapsed time: 2 seconds

Run UMAP/tSNE

GEX.seur <- RunTSNE(GEX.seur, dims=PCs, complexity = 100)
GEX.seur <- RunUMAP(GEX.seur, dims=PCs, n.neighbors = 20, seed.use = 145)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 00:10:26 UMAP embedding parameters a = 0.9922 b = 1.112
## 00:10:26 Read 16562 rows and found 18 numeric columns
## 00:10:26 Using Annoy for neighbor search, n_neighbors = 20
## 00:10:26 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 00:10:28 Writing NN index file to temp file C:\Users\Shaorui\AppData\Local\Temp\RtmpkH5UfV\file4b6c6ed3248
## 00:10:28 Searching Annoy index using 1 thread, search_k = 2000
## 00:10:32 Annoy recall = 100%
## 00:10:32 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 20
## 00:10:33 Initializing from normalized Laplacian + noise (using irlba)
## 00:10:35 Commencing optimization for 200 epochs, with 482664 positive edges
## 00:10:52 Optimization finished
DimPlot(GEX.seur, reduction = "tsne", label = T) + DimPlot(GEX.seur, reduction = "umap", label = T)

DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "preAnno2") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

FeaturePlot(GEX.seur, reduction = "umap", features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rb"))

length(as.vector(unlist(markers.new.s)))
## [1] 92
FeaturePlot(subset(GEX.seur, subset = cnt %in% c("CTL")),
              features = as.vector(unlist(markers.new.s)), cols = c("lightgrey","red"))

FeaturePlot(subset(GEX.seur, subset = cnt %in% c("CKO")),
              features = as.vector(unlist(markers.new.s)), cols = c("lightgrey","red"))

Anno

#
GEX.seur$Anno1 <- as.character(GEX.seur$seurat_clusters)


GEX.seur$Anno1[GEX.seur$Anno1 %in% c(1,2,0,13,12)] <- "EMN1"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(8)] <- "EMN2"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(9,14)] <- "EMN3"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(18,24)] <- "EMN4"

GEX.seur$Anno1[GEX.seur$Anno1 %in% c(3,6,5,4)] <- "IMN1"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(16)] <- "IMN2"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(10)] <- "IMN3"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(23)] <- "IMN4"

GEX.seur$Anno1[GEX.seur$Anno1 %in% c(7)] <- "IN1"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(19)] <- "IN2"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(22)] <- "IN3"

GEX.seur$Anno1[GEX.seur$Anno1 %in% c(15,11)] <- "IPAN1"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(17,20)] <- "IPAN2"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(25)] <- "IPAN3"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(21)] <- "IPAN4"

GEX.seur$Anno1 <- factor(GEX.seur$Anno1,
                         levels = c(paste0("EMN",1:4),
                                    paste0("IMN",1:4),
                                    paste0("IN",1:3),
                                    paste0("IPAN",1:4)))

#
GEX.seur$Anno2 <- as.character(GEX.seur$seurat_clusters)


GEX.seur$Anno2[GEX.seur$Anno2 %in% c(1,2,0,13,12)] <- "EMN1"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(8)] <- "EMN2"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(9,14)] <- "EMN3"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(18)] <- "EMN4"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(24)] <- "EMN5"

GEX.seur$Anno2[GEX.seur$Anno2 %in% c(3,6,5,4)] <- "IMN1"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(16)] <- "IMN2"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(10)] <- "IMN3"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(23)] <- "IMN4"

GEX.seur$Anno2[GEX.seur$Anno2 %in% c(7)] <- "IN1"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(19)] <- "IN2"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(22)] <- "IN3"

GEX.seur$Anno2[GEX.seur$Anno2 %in% c(15)] <- "IPAN1.1"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(11)] <- "IPAN1.2"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(17)] <- "IPAN2.1"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(20)] <- "IPAN2.2"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(25)] <- "IPAN3"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(21)] <- "IPAN4"

GEX.seur$Anno2 <- factor(GEX.seur$Anno2,
                         levels = c(paste0("EMN",1:5),
                                    paste0("IMN",1:4),
                                    paste0("IN",1:3),
                                    paste0("IPAN",c(1.1,1.2,2.1,2.2,3,4))))
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno1") +
  DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno2", label.size = 3.2)

cells.nn <- colnames(GEX.seur)[GEX.seur$preAnno2 %in% c("IPAN3","IPAN4") & GEX.seur$Anno1 == "IMN4"]
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "Anno1", cells.highlight = cells.nn) +
  DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno2", label.size = 3.2)

#

GEX.seur$Anno3 <- as.character(GEX.seur$Anno1)

GEX.seur$Anno3[cells.nn] <- "Mixlike"
GEX.seur$Anno3 <- factor(GEX.seur$Anno3,
                         levels = c(levels(GEX.seur$Anno1),"Mixlike"))
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno3", label.size = 3.2) +
  DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno2", label.size = 3.2)

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "FB.info", cols = color.FB)

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0, group.by = "FB.info", cols = color.FB)

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "Anno3")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0, group.by = "Anno3")

length(cells.nn)
## [1] 42
GEX.seur <- subset(GEX.seur, subset = Anno3 != "Mixlike")
GEX.seur
## An object of class Seurat 
## 24081 features across 16520 samples within 1 assay 
## Active assay: RNA (24081 features, 1046 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno1") +
  DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno2", label.size = 3.2)

DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

FeaturePlot(GEX.seur, reduction = "umap", features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rb"))

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "Anno2")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0, group.by = "Anno2")

markers

# find markers for every cluster compared to all remaining cells, report only the positive ones
Idents(GEX.seur) <- "Anno2"

#GEX.markers.pre <- FindAllMarkers(GEX.seur, only.pos = TRUE, min.pct = 0.05,
#                                  test.use = "MAST",
#                                  #test.use = "wilcox",
#                                  logfc.threshold = 0.25)
GEX.markers.pre <- read.table("Anno2.20240125.markers.csv", header = TRUE, sep = ",")
#GEX.markers.pre %>% group_by(cluster) %>% top_n(n = 8, wt = avg_log2FC)
GEX.markers.pre$cluster <- factor(as.character(GEX.markers.pre$cluster),
                          levels = levels(GEX.seur$Anno2))



markers.pre_t32 <- (GEX.markers.pre %>% group_by(cluster) %>% 
                  filter(pct.1>0.1) %>%
                   top_n(n = 32, wt = avg_log2FC) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene


markers.pre_t48 <- (GEX.markers.pre %>% group_by(cluster) %>% 
                  filter(pct.1>0.05 & gene %in% grep("Rps|Rpl|mt-",GEX.markers.pre$gene,invert = T,value = T)) %>%
                   top_n(n = 48, wt = avg_log2FC) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene

markers.pre_t60 <- (GEX.markers.pre %>% group_by(cluster) %>% 
                  filter(pct.1>0.05 & gene %in% grep("Rps|Rpl|mt-|Gm",GEX.markers.pre$gene,invert = T,value = T)) %>%
                   top_n(n = 60, wt = avg_log2FC) %>%
                    filter(p_val_adj < 0.01) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene

markers.pre_t120 <- (GEX.markers.pre %>% group_by(cluster) %>% 
                  filter(pct.1>0.05 & gene %in% grep("Rps|Rpl|mt-|Gm",GEX.markers.pre$gene,invert = T,value = T)) %>%
                   top_n(n = 120, wt = avg_log2FC) %>%
                    filter(p_val_adj < 0.01) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene
ttt = 557
ttt/60
## [1] 9.283333
ttt/64
## [1] 8.703125
ttt/65
## [1] 8.569231
pp.t60 <- list()

for(i in 1:9){
pp.t60[[i]] <- DotPlot(GEX.seur, features = rev(markers.pre_t60[(64*i-63):(64*i)]))  + coord_flip() + 
           theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
}

pp.t60
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

ttt = 972
ttt/60
## [1] 16.2
ttt/64
## [1] 15.1875
ttt/65
## [1] 14.95385
pp.t120 <- list()

for(i in 1:15){
pp.t120[[i]] <- DotPlot(GEX.seur, features = rev(markers.pre_t120[(65*i-64):(65*i)]))  + coord_flip() + 
           theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
}

#pp.t120
table(GEX.seur$DoubletFinder0.05)
## 
## Singlet 
##   16520
table(GEX.seur$DoubletFinder0.1)
## 
## Singlet 
##   16520
markers.new.ss <- list(EMN=c("Chat","Bnc2","Tox","Ptprt",
                       "Gfra2","Oprk1","Adamtsl1", 
                       "Fbxw15","Fbxw24","Chrna7",
                       "Satb1","Itga6","Cntnap5b",
                       "Pgm5","Chgb",
                       "Nxph1",
                       "Gabrb1","Glp2r","Nebl",
                       "Lrrc7",
                       "Ryr3","Eda",
                       "Hgf","Lama2","Efnb2",
                       "Tac1",
                       "Kctd8",
                       "Ptn",
                       "Ntrk2","Penk","Itgb8",
                       "Fut9","Nfatc1","Egfr","Pdia5",
                       "Ahr","Mgll",
                       "Aff3",
                       "Chrm3"
                       ),
                 IMN=c("Nos1","Kcnab1",
                       "Gfra1","Etv1",
                       "Man1a","Airn",
                       "Adcy2",
                       "Col25a1",
                       "Cmah","Creb5","Vip","Pde1a",
                       "Ebf1","Gpc5"
                       ),
                 IN=c("Npas3","Synpr","St18","Gal",
                      "Nova1",
                      "Cdh10","Kcnk13",
                      "Neurod6","Moxd1","Sctr",
                      "Piezo1","Sst","Adamts9","Kcnn2",
                      "Calb2"),
                 IPAN=c("Calcb","Nmu","Adgrg6","Pcdh10",
                        "Ngfr","Galr1","Il7","Aff2",
                        "Gpr149","Itgb6","Met","Itgbl1",
                        
                        "Cdh6","Cdh8",
                        "Clstn2","Ano2","Ntrk3",
                        "Cpne4","Vwc2l","Cdh9",
                        "Car10","Dcc",
                        "Scgn",
                        "Vcan","Cck","Piezo2","Kcnh7",
                        "Rerg","Bmpr1b","Skap1","Ntng1",
                        "Tafa2","Nxph2"),
                 CONTAM=c("Actl6b","Phox2b"),
                 pDEGs=c("Grid2","Ide","Btaf1","Slc15a2","Ccser1","Hdac9","Rspo2","Grem2"))

pm.All_new.s <- DotPlot(GEX.seur, features = as.vector(unlist(markers.new.ss)), group.by = "Anno2",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All")
pm.All_new.s

pm.CTL_new.s <- DotPlot(subset(GEX.seur,subset=cnt %in% c("CTL")), features = as.vector(unlist(markers.new.ss)), group.by = "Anno2",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="CTL only")
pm.CTL_new.s

pm.CKO_new.s <- DotPlot(subset(GEX.seur,subset=cnt %in% c("CKO")), features = as.vector(unlist(markers.new.ss)), group.by = "Anno2",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="CKO only")
pm.CKO_new.s

cell composition

head(GEX.seur@meta.data)
##                    orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCAAGAAGTCTA-1    CR7d_SI       2107         1263 0.00000000  0.4271476
## AAACCCAAGACATACA-1    CR7d_SI       2345         1327 0.04264392  0.3411514
## AAACCCAAGACTCATC-1    CR7d_SI       5340         2565 0.71161049  0.6179775
## AAACCCAAGAGGTTTA-1    CR7d_SI       2542         1295 0.27537372  0.4327301
## AAACCCAAGGCTCTAT-1    CR7d_SI       2705         1514 1.25693161  0.7024030
## AAACCCAAGGTCTTTG-1    CR7d_SI       1053          765 0.00000000  0.5698006
##                    FB.info      S.Score    G2M.Score Phase RNA_snn_res.2
## AAACCCAAGAAGTCTA-1   CTL.3  0.007075537 -0.006547972     S             6
## AAACCCAAGACATACA-1   CTL.3 -0.002734042  0.015554703   G2M            19
## AAACCCAAGACTCATC-1   CTL.4  0.082419696  0.005693510     S            15
## AAACCCAAGAGGTTTA-1   CTL.1 -0.003888110  0.010933631   G2M            11
## AAACCCAAGGCTCTAT-1   CKO.4  0.008380733  0.021945733   G2M            11
## AAACCCAAGGTCTTTG-1   CKO.3  0.008806639 -0.010166359     S            20
##                    seurat_clusters cnt pANN_0.25_0.005_944 DoubletFinder0.05
## AAACCCAAGAAGTCTA-1               6 CTL           0.0000000           Singlet
## AAACCCAAGACATACA-1              19 CTL           0.1111111           Singlet
## AAACCCAAGACTCATC-1              15 CTL           0.1190476           Singlet
## AAACCCAAGAGGTTTA-1              11 CTL           0.0000000           Singlet
## AAACCCAAGGCTCTAT-1              11 CKO           0.0000000           Singlet
## AAACCCAAGGTCTTTG-1              20 CKO           0.0000000           Singlet
##                    pANN_0.25_0.005_1888 DoubletFinder0.1 sort_clusters preAnno1
## AAACCCAAGAAGTCTA-1            0.0000000          Singlet             7      IMN
## AAACCCAAGACATACA-1            0.0952381          Singlet            20       IN
## AAACCCAAGACTCATC-1            0.1190476          Singlet            15     IPAN
## AAACCCAAGAGGTTTA-1            0.0000000          Singlet            11     IPAN
## AAACCCAAGGCTCTAT-1            0.0000000          Singlet            11     IPAN
## AAACCCAAGGTCTTTG-1            0.0000000          Singlet            21     IPAN
##                    preAnno2 Anno1   Anno2 Anno3
## AAACCCAAGAAGTCTA-1     IMN1  IMN1    IMN1  IMN1
## AAACCCAAGACATACA-1      IN2   IN2     IN2   IN2
## AAACCCAAGACTCATC-1    IPAN1 IPAN1 IPAN1.1 IPAN1
## AAACCCAAGAGGTTTA-1    IPAN1 IPAN1 IPAN1.2 IPAN1
## AAACCCAAGGCTCTAT-1    IPAN1 IPAN1 IPAN1.2 IPAN1
## AAACCCAAGGTCTTTG-1    IPAN2 IPAN2 IPAN2.2 IPAN2
GEX.seur$FB.info <- factor(as.character(GEX.seur$FB.info),
                           levels = c(paste0("CTL.",1:4),
                                      paste0("CKO.",c(1,3,4))))
GEX.seur$rep <- paste0("rep",
                        gsub("CTL|CKO","",as.character(GEX.seur$FB.info)))
DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "FB.info", split.by = "FB.info", 
        ncol = 4, cols = color.FB)

#DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "FB.info", split.by = "FB.info", 
#        ncol = 3, cols = color.FB)
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno1") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno2") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

heatmap.1

cowplot::plot_grid(
pheatmap::pheatmap(table(FB.info=GEX.seur$FB.info,
      clusters=GEX.seur$Anno1),
                   main = "Cell Count",
                   gaps_row = c(4),
      gaps_col = c(4,8,11),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,

pheatmap::pheatmap(100*rowRatio(table(FB.info=GEX.seur$FB.info,
                                clusters=GEX.seur$Anno1)),
                   main = "Cell Ratio",
                   gaps_row = c(4),
      gaps_col = c(4,8,11),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)

barplot.1

stat_Anno1 <- GEX.seur@meta.data[,c("cnt","rep","Anno1")]

stat_Anno1.s <- stat_Anno1 %>%
  group_by(cnt,rep,Anno1) %>%
  summarise(count=n()) %>%
  mutate(prop= count/sum(count)) %>%
  ungroup
## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.
pcom.Anno1 <- stat_Anno1.s %>%
  ggplot(aes(x = Anno1, y = 100*prop, color=cnt)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  theme_classic(base_size = 15) +
  scale_color_manual(values = c("skyblue","pink"), name="") +
  labs(title="Cell Composition of SI data, Anno1", y = "Proportion") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
  
  geom_point(aes(x=Anno1, y = 100*prop, color= cnt),
             position = position_dodge(0.75),
             shape=16,alpha=0.75,size=2.15,
             stroke=0.15, show.legend=F)

pcom.Anno1

sig.test

glm.nb - anova.LRT

Sort.N <- c("CTL","CKO")

glm.nb_anova.Anno1 <- list()

for(x1 in 1:2){
  for(x2 in 1:2){
    if(x2>x1){
      
      stat_Anno1.s_N <- stat_Anno1.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
      
      total.N <- stat_Anno1.s_N %>%
        group_by(cnt, rep) %>%
        summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
      
      rownames(total.N) <- paste0(as.character(total.N$cnt),
                                  "_",
                                  as.character(total.N$rep))
      
      stat_Anno1.s_N$total <- total.N[paste0(stat_Anno1.s_N$cnt,
                                            "_",
                                            stat_Anno1.s_N$rep),"total"]
      
      glm.nb_anova.Anno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
      
      for(AA in levels(stat_Anno1.s_N$Anno1)){
        glm.nb_anova.Anno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
          aaa <- stat_Anno1.s_N %>% filter(Anno1==AA)
          aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
          aaa.anova <- anova(aaa.glmnb, test = "LRT")
          aaa.anova$'Pr(>Chi)'[2]
        }, error=function(e){
        
          tryCatch({
            aaa <- stat_Anno1.s_N %>% filter(Anno1==AA)
            aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
            aaa.anova <- anova(aaa.glmnb, test = "LRT")
            aaa.anova$'Pr(>Chi)'[2]
          }, error=function(e){
            
            tryCatch({
              aaa <- stat_Anno1.s_N %>% filter(Anno1==AA)
              aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
              aaa.anova <- anova(aaa.glmnb, test = "LRT")
              aaa.anova$'Pr(>Chi)'[2]
            }, error = function(e){
              NA
            })     
          })
        })
      }
      glm.nb_anova.Anno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.Anno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]]) 
    }
  }
}
glm.nb_anova.Anno1_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.Anno1)))
rownames(glm.nb_anova.Anno1_df) <- names(glm.nb_anova.Anno1)
colnames(glm.nb_anova.Anno1_df) <- gsub("X","C",colnames(glm.nb_anova.Anno1_df))
glm.nb_anova.Anno1_df
##               EMN1      EMN2      EMN3      EMN4      IMN1      IMN2      IMN3
## CTLvsCKO 0.3344298 0.5904035 0.5503359 0.4107641 0.3560548 0.1619612 0.7252822
##               IMN4          IN1      IN2       IN3      IPAN1     IPAN2
## CTLvsCKO 0.7400615 0.0004540387 0.487195 0.3854356 0.08110179 0.9761208
##              IPAN3    IPAN4
## CTLvsCKO 0.1452842 0.415677
round(glm.nb_anova.Anno1_df,5)
##             EMN1   EMN2    EMN3    EMN4    IMN1    IMN2    IMN3    IMN4     IN1
## CTLvsCKO 0.33443 0.5904 0.55034 0.41076 0.35605 0.16196 0.72528 0.74006 0.00045
##             IN2     IN3  IPAN1   IPAN2   IPAN3   IPAN4
## CTLvsCKO 0.4872 0.38544 0.0811 0.97612 0.14528 0.41568

heatmap.2

cowplot::plot_grid(
pheatmap::pheatmap(table(FB.info=GEX.seur$FB.info,
      clusters=GEX.seur$Anno2),
                   main = "Cell Count",
                   gaps_row = c(4),
      gaps_col = c(5,9,12),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,

pheatmap::pheatmap(100*rowRatio(table(FB.info=GEX.seur$FB.info,
                                clusters=GEX.seur$Anno2)),
                   main = "Cell Ratio",
                   gaps_row = c(4),
      gaps_col = c(5,9,12),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)

barplot.2

stat_Anno2 <- GEX.seur@meta.data[,c("cnt","rep","Anno2")]

stat_Anno2.s <- stat_Anno2 %>%
  group_by(cnt,rep,Anno2) %>%
  summarise(count=n()) %>%
  mutate(prop= count/sum(count)) %>%
  ungroup
## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.
pcom.Anno2 <- stat_Anno2.s %>%
  ggplot(aes(x = Anno2, y = 100*prop, color=cnt)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  theme_classic(base_size = 15) +
  scale_color_manual(values = c("skyblue","pink"), name="") +
  labs(title="Cell Composition of SI data, Anno2", y = "Proportion") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
  
  geom_point(aes(x=Anno2, y = 100*prop, color= cnt),
             position = position_dodge(0.75),
             shape=16,alpha=0.75,size=2.15,
             stroke=0.15, show.legend=F)

pcom.Anno2

sig.test

glm.nb - anova.LRT

Sort.N <- c("CTL","CKO")

glm.nb_anova.Anno2 <- list()

for(x1 in 1:2){
  for(x2 in 1:2){
    if(x2>x1){
      
      stat_Anno2.s_N <- stat_Anno2.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
      
      total.N <- stat_Anno2.s_N %>%
        group_by(cnt, rep) %>%
        summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
      
      rownames(total.N) <- paste0(as.character(total.N$cnt),
                                  "_",
                                  as.character(total.N$rep))
      
      stat_Anno2.s_N$total <- total.N[paste0(stat_Anno2.s_N$cnt,
                                            "_",
                                            stat_Anno2.s_N$rep),"total"]
      
      glm.nb_anova.Anno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
      
      for(AA in levels(stat_Anno2.s_N$Anno2)){
        glm.nb_anova.Anno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
          aaa <- stat_Anno2.s_N %>% filter(Anno2==AA)
          aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
          aaa.anova <- anova(aaa.glmnb, test = "LRT")
          aaa.anova$'Pr(>Chi)'[2]
        }, error=function(e){
        
          tryCatch({
            aaa <- stat_Anno2.s_N %>% filter(Anno2==AA)
            aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
            aaa.anova <- anova(aaa.glmnb, test = "LRT")
            aaa.anova$'Pr(>Chi)'[2]
          }, error=function(e){
            
            tryCatch({
              aaa <- stat_Anno2.s_N %>% filter(Anno2==AA)
              aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
              aaa.anova <- anova(aaa.glmnb, test = "LRT")
              aaa.anova$'Pr(>Chi)'[2]
            }, error = function(e){
              NA
            })     
          })
        })
      }
      glm.nb_anova.Anno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.Anno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]]) 
    }
  }
}
glm.nb_anova.Anno2_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.Anno2)))
rownames(glm.nb_anova.Anno2_df) <- names(glm.nb_anova.Anno2)
colnames(glm.nb_anova.Anno2_df) <- gsub("X","C",colnames(glm.nb_anova.Anno2_df))
glm.nb_anova.Anno2_df
##               EMN1      EMN2      EMN3     EMN4       EMN5      IMN1      IMN2
## CTLvsCKO 0.3344298 0.5904035 0.5503359 0.919722 0.05319255 0.3560548 0.1619612
##               IMN3      IMN4          IN1      IN2       IN3   IPAN1.1
## CTLvsCKO 0.7252822 0.7400615 0.0004540387 0.487195 0.3854356 0.5741269
##             IPAN1.2   IPAN2.1   IPAN2.2     IPAN3    IPAN4
## CTLvsCKO 0.01974064 0.5762993 0.4951306 0.1452842 0.415677
round(glm.nb_anova.Anno2_df,5)
##             EMN1   EMN2    EMN3    EMN4    EMN5    IMN1    IMN2    IMN3    IMN4
## CTLvsCKO 0.33443 0.5904 0.55034 0.91972 0.05319 0.35605 0.16196 0.72528 0.74006
##              IN1    IN2     IN3 IPAN1.1 IPAN1.2 IPAN2.1 IPAN2.2   IPAN3   IPAN4
## CTLvsCKO 0.00045 0.4872 0.38544 0.57413 0.01974  0.5763 0.49513 0.14528 0.41568